The Navstar Global Positioning System, hereafter referred to simply as GPS, is a space based radio positioning network for providing users equipped with suitable receivers highly accurate position, velocity, and time (PVT) information. Developed by the United States Department of Defense (DOD), the space based portion of GPS comprises a constellation of GPS satellites in non-geosynchronous 12 hour orbits around the earth.
Prior art FIG. 1 shows the constellation 100 of GPS satellites 101 in orbit. The GPS satellites 101 are located in six orbital planes 102 with four of the GPS satellites 101 in each plane, plus a number of "on orbit" spare satellites (not shown) for redundancy. The orbital planes 102 of the GPS satellites 101 have an inclination of 55 degrees relative to the equator and an altitude of approximately 20,200 km (10,900 miles) and typically complete an orbit in approximately 12 hours. This positions each of the GPS satellites 101 in such a manner that a minimum of five of the GPS satellites 101 are normally observable (above the horizon) by a user anywhere on earth at any given time.
GPS position determination is based upon a concept referred to as time of arrival (TOA) ranging. Each of the orbiting GPS satellites 101 broadcasts spread spectrum microwave signals encoded with positioning data and satellite ephemeris information. The signals are broadcast on two frequencies, L1 at 1575.42 MHz and L2 at 1227.60 MHz, modulated using bi-phase shift keying techniques. Essentially, the signals are broadcast at precisely known times and at precisely known intervals. The signals are encoded with their precise time of transmission. A user receives the signals with a GPS receiver designed to time the signals and to demodulate the satellite orbital data contained in the signals. Using the orbital data, the GPS receiver determines the time between transmission of the signal by the satellite and reception by the receiver. Multiplying this by the speed of light gives what is termed the pseudo range measurement of that satellite. If the GPS receiver clock were perfect, this would be the range measurement for that satellite, but the imperfection of the clock causes it to differ by the time offset between actual time and receiver time. Thus, the measurement is called a pseudo range, rather than a range. However, the time offset is common to the pseudo range measurements of all the satellites. By determining the pseudo ranges of four or more satellites, the GPS receiver is able to determine its location in three dimensions, as well the time offset. Thus, a user equipped with a proper GPS receiver is able to determine his PVT with great accuracy, and use this information to navigate safely and accurately from point to point, among other uses.
As described above, the broadcast signals from the GPS satellites include satellite ephemeris information, essentially a set of orbital parameters, in addition to a precise time of transmission. The ephemeris information from each respective satellite provides a precise description of that satellite's orbit, referred to as the ephemeris, and more simplified descriptions of all the GPS satellite orbits, referred to as almanacs. Each ephemeris precisely characterizes its respective satellite's path by giving the parameters of a Keplerian orbit and several associated corrections to that orbit. This information is required to determine accurate pseudo range measurements to the respective satellite. Generally, the ephemeris information is used to calculate the positions of the GPS satellite with respect to earth-centered earth-fixed coordinates, and in turn, with respect to the GPS receiver. These positions are used in the TOA ranging calculations to determine the pseudo ranges and the time offset.
Specifically, the ephemeris for a given satellite comprises:
a) The square-root of the semi-major orbital axis, which is half the greatest dimension of the elliptical orbit. The square-root ranges from 0 to about 8192 square-root-meters, corresponding to a semi-major axis a of about 67 megameters, with a resolution corresponding to about 20 millimeters for the typical orbit that has a semi-major axis of 26560 kilometers. Transmission of the square-root rather than the semi-major axis itself allows calculation of the mean motion n.sub.0, which is the average orbital angular velocity, without using a square-root operation: EQU n.sub.0 =(.mu./a.sup.3).sup.0.5 =.mu..sup.0.5 /(a.sup.0.5).sup.3
where .mu. is the earth's gravitational parameter, commonly the WGS-84 value of 398600.5 kilometers-cubed-per-second-squared. The value of n.sub.0 is typically 146 microradians-per-second, corresponding to about 3874 meters-per-second.
b) The ephemeris reference time of the week T.sub.0, which is the time from the beginning of the current GPS week to a time near the middle of the time interval for which the ephemeris is valid. The reference time T.sub.0 is an integer multiple of 16 seconds. Time T.sub.0 is commonly subtracted from the given time of interest T to get ephemeris time t, that is, time referred to the ephemeris's epoch: EQU t=T-T.sub.0
c) The mean-motion difference .increment.n, which is a correction to the mean motion no as computed above: EQU n=n.sub.0 .increment.n
The range of .increment.n corresponds to about .+-.311 millimeters-per-second with a resolution of about 9.5 micrometers-per-second for a typical orbit.
d) The reference-time mean anomaly M.sub.0, which is the angle of the satellite at the reference time T.sub.0 from perigee which would result from uniform angular motion. Perigee is the point in the orbit closest to the earth. The resolution of M.sub.0 corresponds to about 39 millimeters for a typical orbit. M.sub.0 is used to find the satellite's mean anomaly M at ephemeris time t: EQU M=M.sub.0 nt
e) The orbital eccentricity e, which is the fraction of the semi-major axis by which the earth is distant from the orbital center. The dimensionless fraction e ranges from 0 to 0.03 with a resolution corresponding to about 3 millimeters for a typical orbit. Kepler's equation relates mean anomaly M to eccentric anomaly E, which is the actual at the reference time T.sub.0 about the orbital center from perigee to the satellite: EQU M=E-e sin E
Kepler's equation is commonly solved iteratively to get eccentric anomaly E from mean anomaly M. True anomaly v, which is the actual angle about the earth's center from perigee to the satellite, is commonly calculated from eccentric anomaly E by any of various formulas, for example: ##EQU1##
f) The argument of perigee .omega., which is the angle about the earth's center from the ascending node to perigee. The ascending node is the intersection of the orbit with the plane of the equator at which the satellite passes from the southern to the northern hemisphere. The resolution of .omega. corresponds to about 39 millimeters for a typical orbit. The argument of perigee .omega. is commonly added to the true anomaly v to get the uncorrected argument of latitude .PHI., which is the angle about the earth's center from the ascending node to the satellite: EQU .PHI.=.omega.+v
g) The second-harmonic correction coefficients C.sub.uc and C.sub.us to the argument of latitude .PHI.. The ranges of C.sub.uc and C.sub.us correspond to about .+-.1621 meters with a resolution of about 49 millimeters for a typical orbit. They multiply the cosine and sine respectively of twice the uncorrected argument of latitude .PHI. to correct the argument of latitude u: EQU u=.PHI.+C.sub.us cos 2.PHI.+C.sub.us sin 2.PHI.
h) The second-harmonic correction coefficients C.sub.rc and C.sub.rs to the radius r. The ranges of C.sub.rc and C.sub.rs are about .+-.1024 meters with a resolution of about 31 millimeters. They multiply the cosine and sine respectively of twice the uncorrected argument of latitude .PHI. to correct the radius r: EQU r=a(1-e cos E)+C.sub.rc cos 2.PHI.+C.sub.rs sin 2.PHI.
The argument of latitude u and the radius r are commonly combined to calculate the satellite's position in the orbital plane (x',y'): EQU x'=r cos u EQU y'=r sin u
i) The orbital inclination i.sub.0 at the reference time, which is the angle that the orbital plane makes with the plane of the equator, typically about 55 degrees. The resolution of i.sub.0 corresponds to about 39 millimeters for a typical orbit.
j) The rate of orbital inclination i. The range of i corresponds to .+-.155 millimeters-per-second with a resolution of about 9.5 micrometers-per-second for a typical orbit.
k) The second-harmonic correction coefficients C.sub.ic and C.sub.is to the inclination i. The ranges of C.sub.ic and C.sub.is correspond to about .+-.1621 meters with a resolution of about 49 millimeters for a typical orbit. They multiply the cosine and sine respectively of twice the uncorrected argument of latitude .PHI. to correct the inclination i: EQU i=i.sub.0 +it+C.sub.uc cos 2.PHI.+C.sub.us sin 2.PHI.
l) The right ascension of the ascending node at the reference time .OMEGA..sub.0, which is the angle measured eastward from the vernal equinox to the ascending node. The resolution of .OMEGA..sub.0 corresponds to about 39 millimeters for a typical orbit.
l) The rate of right ascension .OMEGA.. The range corresponds to about .+-.80 meters-per-second with a resolution of about 9.5 micrometers-per-second for a typical orbit. .OMEGA..sub.0 and .OMEGA. are commonly combined with .OMEGA..sub.e, the WGS-84 earth's rotation rate, to calculate the ascending node's longitude: EQU .OMEGA.=.OMEGA..sub.0 .OMEGA.t-.OMEGA..sub.e T
The ascending node's longitude .OMEGA. and the inclination i are commonly combined with the satellite's position in the orbital plane to get the satellite's position in earth-centered earth-fixed coordinates (x, y, z): EQU x=x' cos .OMEGA.-y' cos i cos .OMEGA. EQU y=x' sin .OMEGA.+y' cos i cos .OMEGA. EQU z=y' sin i
The above calculations describe the GPS satellite orbit as a function of time. In other words, for a given time, the equations yield the position of a given GPS satellite in its orbit. As described above, the positions of the various GPS satellites are required in order to compute pseudo ranges. Time is used as an input variable to the above equations. Using the above equations and parameters transmitted from the GPS satellites, a computer system (e.g., embedded within a GPS receiver) computes time referenced position fixes for the GPS satellites. The above equations are standardized in a GPS technical specification, "Navstar GPS Space Segment/Navigation User Interface," ICD-GPS-200 Revision A, Table 20-IV, Rockwell International Corporation, Downey, Calif., 1984, and are well known and widely used. Those desiring additional or more detailed information are referred to the above document, which is incorporated herein as background material.
The problem with the above ICD-GPS-200 equations is that the iterative solution of Kepler's equation (upon which the ICD-GPS-200 equations are based) and the many transcendental functions involved require more computation than is practical for many applications. For example, a typical GPS navigation receiver in an aircraft needs to compute an updated GPS position of the aircraft at least, on average, once per second. This requires the computation of pseudo ranges from at least four GPS satellites, which in turn, requires the determination of the positions of those GPS satellites at least once per second. In more demanding cases, for example, where the GPS position of the satellite is updated several times per second, the computational load of iteratively solving for the positions of the GPS satellites using the ICD-GPS-200 equations may be unworkable.
One technique for reducing this computational load is interpolation. The interpolation technique involves solving for samples of GPS satellite position over a relatively long period of time, and then determining interpolated GPS satellite positions between the solved positions. Interpolation involves a trade-off among position accuracy, frequency of calculation, and the complexity of the interpolation, which is commonly in the form of a polynomial in time. Generally, for a given GPS satellite, several samples are calculated using the ICD-GPS-200 equations. A polynomial curve is then numerically fitted to the samples. This curve allows positions between the samples to be determined through interpolation. The GPS satellite's position at a particular time is determined using the polynomial. In so doing, numerous satellite positions can be determined even though only a much smaller number of samples are actually computed. This greatly reduces the computational load required for maintaining high update rates in GPS navigation systems. Accordingly, the interpolation technique is well established within the GPS industry and widely used. Most GPS navigation systems implement the interpolation technique for providing GPS position fixes.
Another technique for reducing the computational load is numerical integration. Numerical integration avoids the transcendental functions and iteration completely, except for an initialization. Numerical integration also involves a trade-off, among position accuracy, frequency of iteration, and the complexity of the integration step. In order to hold reasonable accuracy over the one-hour duration of even the shortest GPS ephemeris, conventional integration must iterate far too often to be at all practical. Accordingly, the use of the numerical integration is not currently supported in current GPS navigation systems.
It should be appreciated that although the numerical integration technique is not currently used for computing GPS satellite samples, an unconventional form of numerical integration, "Leapfrog" integration, has advantages for a lossless system, like a Keplerian orbit. For examples and more detailed discussion, refer to Robert Leonard Nelson, Jr., and Hercule Kwan, "`Leapfrog` Methods for Numerical Solution of Differential Equations, Sinewave Generation, and Magnitude Approximation," IEEE Asilomar Conference, October, 1995. In general, Leapfrog numerical integration updates position and velocity on alternate half-intervals rather than aligned in time: EQU Increment p by f(p).delta.t=-.mu.p.parallel.p.parallel..sup.-3 .delta.t EQU Increment p by p .delta.t
where p is the position vector at the ends of the intervals, p is velocity at the midpoints, .mu. is the earth's gravitational constant, and .delta.t is the time step.
Leapfrog integration accumulates no radial error for circular motion, as is the case with a perfectly circular orbit. Furthermore, a simple update-multiplier adjustment eliminates along-orbit error: EQU Decrement p by 2(p.OMEGA.) sin (.OMEGA. .delta.t/2) EQU Increment p by 2(p/.OMEGA.) sin (.OMEGA. .delta.t/2)
where .OMEGA. is angular velocity and .OMEGA. .delta.t/2 is the angle step.
A re-formulation avoids multiplication: EQU Decrement (p/.OMEGA.) by 2p sin (.OMEGA. .delta.t/2) EQU Increment p by 2(p/.OMEGA.) sin (.OMEGA. .delta.t/2)
or
______________________________________ Decrement .sub.p.sup.+ by p/d Increment p by .sub.p.sup.30 /d ______________________________________
where ##STR1## is the "pseudo-velocity", proportional to velocity but equal dimensionally and in magnitude to the radius, and d is an integer power of two which now determines .delta.t: EQU .delta.t=(2/.OMEGA.) arc sin (1/2d)
Although leapfrog integration can be made exact for circular motion or for any motion that is a sinusoidal function of the independent variable, the GPS satellite's position in its orbit is not generally exactly circular nor a sinusoidal function of time. Thus, numerical integration is subject to rapidly accumulating error as the step size increases. To maintain acceptable accuracy, the step size needs to be relatively small, which incurs the corresponding computer load penalty. For these reasons, numerical integration is not practical for GPS navigation systems.
Thus, what is required is a highly accurate method for calculating the positions of orbiting GPS satellites. What is further required is a computationally efficient method of calculating orbital position which provides high accuracy with a minimal demand for computation resources. The required method should provide for the computation of satellite positions in less time, or with a less powerful computer system, or both. What is further required is a method for calculating the positions of orbiting GPS satellites, or any satellite, which does not lose accuracy for larger step sizes, allowing the step size to be as large as the sample spacing wanted for interpolation. The present invention provides a novel solution to the above requirements.